Przed oddaniem zadania upewnij się, że wszystko działa poprawnie. Uruchom ponownie kernel (z paska menu: Kernel$\rightarrow$Restart) a następnie wykonaj wszystkie komórki (z paska menu: Cell$\rightarrow$Run All).

Upewnij się, że wypełniłeś wszystkie pola TU WPISZ KOD lub TU WPISZ ODPOWIEDŹ, oraz że podałeś swoje imię i nazwisko poniżej:

---

Zmienne kategoryczne w modelach liniowych¶

Podejść do "wektoryzacji" zmiennych kategorycznych jest wiele i ich stosowanie jest niezbędne w przypadku większości modeli. W przypadku modeli liniowych dobór odpowiedniego podejścia jest jednak wyjątkowo istotny, ponieważ model nie odtworzy ukrytych zależności - w przeciwieństwie do np. sieci neuronowych (przynajmniej w teorii). Oto wybór możliwości:

1. Osobny model dla każdej z kategorii

   Mało wydajne, ponieważ zwiększamy liczbę stopni swobody (VC-dimensiom, klątwa wymiarowości itp) przy tej samej liczbie danych. Ale jeżeli każda z klas posiada inny proces generujący, to po co łączyć w jeden model?

2. Utworzenie $k$ bądź $k - 1$ nowych zmiennych niezależnych (gdzie $k$ to liczba kategorii)

   Popularne i łatwe w użyciu podejście, nie wymaga zrozumienia zbioru danych, daje dobre wyniki dla modeli nie-liniowych. Można przeprowadzić na wiele sposobów - https://stats.idre.ucla.edu/spss/faq/coding-systems-for-categorical-variables-in-regression-analysis-2/

3. Warunkowanie istniejących zmiennych niezależnych

   Wymaga zrozumienia zbioru danych (dla jakich zmiennych niezależnych wpływ na zmienną zależną może być różny zależnie od kategorii), ale dla modeli liniowych będzie dawać najlepsze wyniki (jeżeli przeprowadzone prawidłowo = zgodnie z procesem generatywnym/z jego najlepszą liniową aproksymacją).

   Poza trudnością w identyfikacji warunkowania, pojawiają się również trudności implementacyjne. Ręczne przygotowanie zbioru danych w ten sposób jest czasochłonne i łatwo o błąd. Bardzo dobrym rozwiązaniem jest podejście z pakietu R, tzw. formula function, którą posiada bardzo funkcjonalną składnię. W Pythonie z tej samej składni można skorzystać w pakiecie statsmodels.

Zastosujemy rozwiązanie nr 2. w wersji dummy / one-hot / one-of-K coding, czyli każda z klas będziała miała swoją kolumnę wypełnioną zerami i jedynkami. Dla jednej z klas (klasy referencyjnej) moglibyśmy nie tworzyć kolumny, ponieważ wyraz wolny "przejąłby" referencyjny wpływ na zmienną zależną. Gdybyśmy stosowali model inny niż liniowy z wyrazem wolnym takie podejście mogłoby się nie sprawdzić. Wyraz wolny, wkrtóce zniknie z naszego modelu liniowego, więc pozostawiamy $k$ kolumn.

A teraz, wcześnie dopasowany na danych treningowych transformator należy zaaplikować na danych testowych:

Oryginalną kolumnę zostawiamy w DataFrameach, ponieważ ułatwi nam to wizualizację.

Ridge Regression¶

Model: $$ y_i = \pmb{x}_i \pmb{\beta} $$

Regularyzacja $L_2$:

$$ \mathcal{R}(\pmb{\beta}) = \frac{\lambda}{2} \left\|\pmb{\beta}\right\|_2^2,$$

gdzie $\lambda$ to parametr metody regulujący siłę kary za złożoność.

Funkcja kosztu: $$ S(\pmb{\beta}) = \frac{1}{2} \sum_{i=1}^N \left| y_i - (\beta_0 + \pmb{x}_i \pmb{\beta}) \right|^2 + \frac{\lambda}{2} \left\|\pmb{\beta}\right\|^2, $$ gdzie macierz $X$ nie posiada tutaj dodatkowej kolumny $1$ odpowiadającej za wyraz wolny. Jest on traktowany osobno, ponieważ nie chcemy na niego zakładać regularyzacji.

Estymacja:

1. Analityczna Najpierw musimy dopasować wyraz wolny $\beta_0$. Ma być to wartość o najmniejszym błędzie średniokwadratowym, gdyby wszystkie cechy miały wartość 0.

   $$ \hat{\beta}_0 = \frac{1}{N} \sum_{i=1}^N y_i $$

   Następnie wyraz wolny należy odjąć od zmiennej niezależnej co daje następującą funkcję kosztu: $$ S(\pmb{\beta}) = \frac{1}{2} \sum_{i=1}^N \left| (y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta}) \right|^2 + \frac{\lambda}{2} \left\|\pmb{\beta}\right\|^2, $$ którą można różniczkować względem parametrów $\pmb{\beta}$ i przyrównać do $0$. Dostajemy następujący estymator $$ \hat{\pmb{\beta}} = (\lambda I_N + X^\intercal X)^{-1}X^\intercal (Y - \hat{\beta}_0) $$

2. Iteracyjna - Gradient Descent

   W tym podejściu do $S(\pmb{\beta})$ podchodzi się jak do funkcji kosztu nie wchodząc w jej interpretację i własności. Jest to problem z zakresu Convex Optimization.

Standaryzacja zmiennych¶

Standaryzowania zmiennych dokonuje się poprzez odjęcie średniej i przeskalowanie do jednostkowej wariancji.

Dla próbki $x$ ustandaryzowana wartość $z$ obliczana jest w następujący sposób:

$$ z = \frac{x - u}{s},$$

gdzie $u$ i $s$ to odpowiednio średnia i odchylenie standardowe populacji.

Dobrą praktyką jest standaryzować (ewentualnie skalować do przedziału) zmienne niezależne oraz zależne przed stosowaniem modeli. Dotychczas tego nie robiliśmy, ponieważ

• dla modelu regresji liniowej skala zmiennych nie ma znaczenia
• mieliśmy jeden model.

W poprzednim zeszycie wspomniane zostało, że przy regularyzowanej regresji skalowanie jest niezbędne - skoro nakładamy taką samą karę na wszystkie parametry, to powinny one mieć wpływ w tej same jednostce na wszystkie zmienne niezależne.

Ponadto dla ustandaryzowanych zmiennych łatwiej jest porównywać modele między sobą, a także dobierać rozkłady prior dla parametrów, a tym będziemy się zajmować w tym zeszycie.

Należy jednak pamiętać, że współczynniki dla ustandaryzowanych zmiennych są cięższe do interpretacji. Dlatego wiele pakietów statystycznych domyślnie dokonuje standaryzacji i raportuje parametry dla ustandaryzowanych zmiennych, a także przeskalowane do oryginalnych rozkładów cech.

W tym zeszycie intercept nie będzie traktowany, jak kolejna z cech. Nie chcemy go także standaryzować - otrzymalibyśmy wektor zer.

Pewien niepokój powinno budzi standaryzowanie zmiennych powstałych w ramach one-hotowania. Pomimo, że mają wartości liczbowe, to ciężko o nich myśleć jako o zmiennych ciągłych. Z drugiej strony chcemy zachować tą samą skalę dla wszystkich zmiennych w celu regularyzacji. Nie istnieje jednoznacza odpowiedź na to pytanie. Większość debaty w Internecie sprowadza się do

sklearn zwróci błąd jeżeli wywołamy StandardScaler przed OneHotEncoder

Należy pamiętać, że jesteśmy teraz poza światem założeń, więc najlepiej zewaluować oba podejścia w różnych warunkach i na tej podstawie wyciągnąć wnioski.

Zastosujemy implementację Ridge Regression z pakietu sklearn, która estymuje wyraz wolny w pkt 1. powyżej, czyli nie regularyzuje go. Posiada ona następujące solvery (dokumentacja):

• ‘auto’ chooses the solver automatically based on the type of data. ‘svd’ uses a Singular Value Decomposition of X to compute the Ridge coefficients. More stable for singular matrices than ‘cholesky’.

• ‘cholesky’ uses the standard scipy.linalg.solve function to obtain a closed-form solution.

• ‘sparse_cg’ uses the conjugate gradient solver as found in scipy.sparse.linalg.cg. As an iterative algorithm, this solver is more appropriate than ‘cholesky’ for large-scale data (possibility to set tol and max_iter).

• ‘lsqr’ uses the dedicated regularized least-squares routine scipy.sparse.linalg.lsqr. It is the fastest and uses an iterative procedure.

• ‘sag’ uses a Stochastic Average Gradient descent, and ‘saga’ uses its improved, unbiased version named SAGA. Both methods also use an iterative procedure, and are often faster than other solvers when both n_samples and n_features are large. Note that ‘sag’ and ‘saga’ fast convergence is only guaranteed on features with approximately the same scale. You can preprocess the data with a scaler from sklearn.preprocessing.

Wśród nich znajdziemy "analityczne" rozwiązanie w oparciu o operacje algebraiczne, a także interacyjne.

$\lambda$ nazywa się tutaj alpha :)

Od regularyzacji do Bayesian Linear Regression¶

(Zakładamy tutaj, że mamy już estymator wyrazu wolnego $\hat{\beta}_0$)

Przypomnijmy, że regresja liniowa to $$ (y_i - \hat{\beta}_0) = \pmb{x}_i \pmb{\beta} + \epsilon_i, $$ gdzie zakładamy $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$. Czyli równoważnie możemy zapisać ten model jako: $$ (y_i - \hat{\beta}_0) | \pmb{x_i} \sim \mathcal{N}(\pmb{x}_i \pmb{\beta}, \sigma^2). $$

Log-likelihood takiego modelu to:

$$ \log p(D|\pmb{\beta}) = \frac{1}{N} \sum_{i=1}^N \log \left[ \frac{1}{\sqrt{2 \pi}\sigma}\exp \left( - \frac{((y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta})^2}{2 \sigma^2} \right)\right] = \text{const} - \frac{1}{2N\sigma^2}\sum_{i=1}^N \left| (y_i - \hat{\beta}_0) - \pmb{x}_i \pmb{\beta}) \right|^2 $$

Gdy spojrzymy na Ridge Regression od strony MAP (nawiązanie do ostatnich zajęć) $$ \arg\max_{\pmb{\beta}} \log p(\pmb{\beta}| D) = \arg\max_{\pmb{\beta}} [ \log p(\pmb{\beta}) + \log p (D | \pmb{\beta})], $$ to okazuje się, że Ridge Regression to probabilistyczny model ze standardowym Gaussowskim priorem na wagach.

Niech $\pmb{\beta} \sim \mathcal{N}(\pmb{m}, V)$: $$ \log p(\pmb{\beta}) = - \frac{1}{2} (\pmb{\beta} - \pmb{m})^\intercal V^{-1} (\pmb{\beta} - \pmb{m}) + \text{constant},$$ gdy $\pmb{m}=0$ i $V = I_{m}$ to $$ \log p(\pmb{\beta}) = - \frac{1}{2} \|\pmb{\beta}\|_2^2 + \text{constant},$$ czyli otrzymaliśmy Ridge Regression!

Bayesian Linear Regression¶

Mamy zatem taki sam model, jak w przypadku Linear Regression, ale dla parametrów wybieramy rozkład prior. Na początek załóżmy, że wariancja błędu jest znana i skupmy się współczynnikach. Likelihood modelu wyrażony jest w następujący sposób $$ p(Y|X, \pmb{\beta}, \beta_0, \sigma^2) = \mathcal{N}(Y | \beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N) \propto \exp \left(- \frac{1}{2\sigma^2} (Y - \beta_0 \pmb{1}_N - X \pmb{\beta})^\intercal (Y - \beta_0 \pmb{1}_N - X \pmb{\beta}) \right) ,\tag{1}$$ gdzie przez $\mathcal{N}(Y | \beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N)$ oznaczamy prawdopobieństwem obserwacji $Y$ przy rozkładzie normalnym o parametrach $(\beta_0 + X \pmb{\beta}, \sigma^2 \mathrm{I}_N)$.

Na parametr $\beta_0$ (wyraz wolny) zakładamy nieinformatywny prior, czyli jego estymacja sprowadza się do MLE i tak jak wcześniej $$ \hat{\beta}_0 = \frac{1}{N} \sum_{i=1}^N y_i = \bar{Y}.$$ Dla ułatwienia notacji będziemy zakładać, że $Y$ zostało wycentrowane (wtedy możemy pominąć wyraz wolny), zatem $$ Y \leftarrow Y - \bar{Y}.$$

Wyśrodkujmy zmienną niezależną w naszym zbiorze danych.

Zatem $(1)$ upraszacza się do $$ p(Y|X, \pmb{\beta}, \sigma^2) \propto \exp \left(- \frac{1}{2\sigma^2} (Y - X \pmb{\beta})^\intercal (Y - X \pmb{\beta}) \right)$$ które w dalszym ciągu jest Gaussianem, więc rozkład sprzężony to również Gaussian, oznaczmy go w następujący sposób $$ p(\pmb{\beta}) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_0, \pmb{V}_0)$$

Wtedy posterior obliczamy jako $$ p(\pmb{\beta}| X, Y, \sigma^2) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_0, \pmb{V}_0) \mathcal{N}(Y | X\pmb{\beta}_0, \sigma^2 \mathrm{I}_N) = \mathcal{N}(\pmb{\beta} | \pmb{\beta}_N, \pmb{V}_N),\tag{2}$$ gdzie $$ \pmb{\beta}_N = \pmb{V}_N \pmb{V}_0^{-1} \pmb{\beta}_0+\frac{1}{\sigma^2}\pmb{V}_N X^\intercal Y $$ $$ \pmb{V}_N = \sigma^2 (\sigma^2 \pmb{V}_0^{-1} + X^\intercal X)^{-1} $$ $$ \pmb{V}_N^{-1} = \pmb{V}_0^{-1} + \frac{1}{\sigma^2} X^\intercal X $$

Skąd to się wzięło? Murphy:eq4.125, Bishop:eq2.116

• Jeżeli $X$ jest ustandaryzowany, to rozsądnym jest $\pmb{\beta}_0 = \pmb{0}$
• Jeżeli na tym etapie chcemy przestać być "Bayesowscy", to możemy zrobić MAP - tutaj znamy pełen rozkład (bo Gaussian), więc bierzemy modę (która równa jest średniej dla tego rozkładu)

Ale my nie chcemy opuszczać "Bayesowskiej" drogi i chcemy doprowadzić do znalezienia rozkładu predykcyjnego (ang. predictive distribution). Możemy tak zrobić, ponieważ znamy pełen posterior, a nie tylko funkcję do niego wprost proporcjonalną.

Przez $\mathcal{D}=(X, Y)$ oznaczmy dotychczasowy zbiór treningowy, składający się z $N$ par. Chcemy znaleźć rozkład (predykcyjny) $y$ dla nowej próbki $\pmb{x}$. W ogólność chcemy wykonać następującą operację $$ p(y|\pmb{x}, \mathcal{D}) = \int p(y|\pmb{x}, \pmb{\beta}) p(\pmb{\beta}|\mathcal{D}) \mathrm{d} \pmb{\beta},$$ gdzie $p(\pmb{\beta}|\mathcal{D})$ to posterior (policzony na bazie priora i likelihoodu). W naszym wypadku jest dodatkowo sparametryzowany przez $\sigma^2$.

Zatem finalnie

$$ p(y|\pmb{x}, \mathcal{D}, \sigma^2) = \int \mathcal{N}(y | \pmb{x}\pmb{\beta}, \sigma^2) \mathcal{N}(\pmb{\beta} | \pmb{\beta}_N, \pmb{V}_N) \mathrm{d} \pmb{\beta} = \mathcal{N}\left(y | x\pmb{\beta}_N, \sigma^2_N(\pmb{x})\right),\tag{3}$$

gdzie $$\sigma^2_N(\pmb{x}) = \sigma^2 + \pmb{x}\pmb{V}_N\pmb{x} $$

Warto zwrócić uwagę na $\sigma^2_N(\pmb{x})$, w której założona wariancja błędu jest powiększona o dodatkowy czynnik, który dodaje wariancji od posteriora wektora wag.

W przypadku probabilistycznego modelu nie mówimy już o przedziałach ufności i przedziałach predykcyjnych, tylko o credible intervals wynikających wprost z rozkładu predykcyjnego.

UWAGA 1: Taki model jest idealny dla uczenia przyrostowego. Dobry opis Bishop:p.154-156

UWAGA 2: Możemy znaleźć posterior dla modelu z priorem na $\sigma^2$ -> Murphy:p.236 rozdział 7.6.3 Bayesian inference when $\sigma^2$ is unknown

Empirical Bayes / Evidence Procedure¶

W powyższych rozważaniach mamy dwa hiper-parametry:

• $\sigma^2$ - wariancja błędu; częściej mówi się o parametrze precision $\lambda=1/\sigma^2$
• $\pmb{V}_0$ - macierz kowariancji dla priora na współczynnikach; najczęściej zakłamy sferyczny prior, więc $\pmb{V}_0 = \alpha^{-1} \mathrm{I}_N$ skupimy się na kolejnym precision $\alpha$

Oznaczmy $\eta= (\alpha, \lambda)$.

Popularnym rozwiązaniem tego problemu byłby Randomized Grid K-CV Search, gdzie dla losowego podzbioru konfiguracji z siatki przeprowadzona byłaby K-cross walidacja.

Alternatywą jest EB. Aby być w pełni Bayesowscy powinniśmy nałożyć hyperpriors na $\eta$. Wtedy możemy obliczyć posterior $\eta$ i rozkład predykcyjny ma następującą postać: $$ p(y|\pmb{x}, \mathcal{D}) = \int p(y|\pmb{x}, \pmb{\beta}) p(\pmb{\beta}|\mathcal{D}, \alpha, \lambda) p(\alpha, \lambda|\mathcal{D}) \mathrm{d} \pmb{\beta} \mathrm{d} \alpha \mathrm{d} \lambda$$

(Hiper-priory też mogą mieć swoje parametry i tak dalej...) Wycałkowanie $\alpha, \lambda$ jest analitycznie niemożliwe, obliczeniowo również. Pierwsze rozwiązanie to MAP - jeżeli posterior jest wyraźnie skupiony wokół wartości $\hat{\eta}$, to wybieramy ją i stosujemy jako plug-in estimate. By zaproponować alternatywę przypomnijmy jaką postać ma posterior $$ p(\alpha, \lambda | \mathcal{D}) \propto p(\mathcal{D}|\alpha, \lambda) p(\alpha, \lambda) $$ Jeżeli prior jest "płaski" i nie daje posteriora o skupionej masie, możemy zwrócić się ku MLE na likelihoodzie $p(\mathcal{D}|\alpha, \lambda)$. W obu przypadkach należy obliczyć posterior $$p(\mathcal{D}|\alpha, \lambda) = \int p(\mathcal{D}|\pmb{\beta}, \lambda) p(\pmb{\beta}|\alpha) \mathrm{d} \pmb{\beta},$$ a w przypadku MLE także przeprowadzić jego maksymalizację (opiera się na Expectation Maximization, które omówione zostanie w dalszej części semestru). Obie te rzeczy są dobrze opisane w Bishop:p.166-169 (UWAGA - inna notacja), rozdziały:

• 3.5.1 Evaluation of the evidence function
• 3.5.2 Maximizing the evidence function

UWAGA: W tym przypadku online learning nie jest wprost aplikowalny - "empiryczna część" potrzebuje próbki by dokonać maksymalizacji (posteriora bądź evidence). Można spróbować Expectation Maximization w konfiguracji batchowej - ciekawy temat na projekt.

Implementacja z rozkładem sprzężonym na $\pmb{\beta}$ i znanym $\alpha$ oraz $\lambda$¶

Zadanie 1 (2pkt)¶

W klasie poniżej zaimplementuj metodę predict tak by zwracała rozkład predykcyjny $Y$ dla zadanych danych wejściowych x. Wszystkie potrzebne wzory to $(2)$ i $(3)$, natomiast wszystkie potrzebne, nieznane parametry modelu zostały przypisane do zmiennych w konstruktorze oraz metodzie fit.

UWAGA: dodaj wyraz wolny do średniej

Implementacja z rozkładem sprzężonym na $\pmb{\beta}$ i nieznanym $\alpha$ oraz $\lambda$¶

Jak wspomnieliśmy powyżej ta metoda może opierać się na:

1. MLE na likelihoodzie (Evidence Maximization) Wyprowadzenie (na podstawie Bishop) oraz implementację tego podejścia można znaleźć pod tym linkiem - http://krasserm.github.io/2019/02/23/bayesian-linear-regression/

2. MAP Rozkładem sprzężonym dla precision jest rozkład Gamma, który zależnie od parametrów preferuje małe bądź większe wartości. Implementacja tego podejścia jest dostępna w pakiecie sklearn pod nazwą sklearn.linear_model.BayesianRidge. Wykorzystamy ją tutaj wraz z omówieniem parametryzacji.

   UWAGA: tutaj $\lambda$ i $\alpha$ są zamienione!

Wybrane parametry klasy BayesianRidge:

• n_iter: int, default=300

  Maximum number of iterations. Should be greater than or equal to 1.

• tol: float, default=1e-3

  Stop the algorithm if w has converged.

• alpha_1: float, default=1e-6

  Hyper-parameter : shape parameter for the Gamma distribution prior over the alpha parameter.

• alpha_2: float, default=1e-6

  Hyper-parameter : inverse scale parameter (rate parameter) for the Gamma distribution prior over the alpha parameter.

• lambda_1: float, default=1e-6

  Hyper-parameter : shape parameter for the Gamma distribution prior over the lambda parameter.

• lambda_2: float, default=1e-6

  Hyper-parameter : inverse scale parameter (rate parameter) for the Gamma distribution prior over the lambda parameter.

• `alpha_init: float, default=None

  Initial value for alpha (precision of the noise). If not set, alpha_init is 1/Var(y).

     New in version 0.22.

• lambda_init: float, default=None

  Initial value for lambda (precision of the weights). If not set, lambda_init is 1.

     New in version 0.22.

• fit_intercept: bool, default=True

  Whether to calculate the intercept for this model. The intercept is not treated as a probabilistic parameter and thus has no associated variance. If set to False, no intercept will be used in calculations (i.e. data is expected to be centered).

Ta implementacja nie wspiera inkrementalnego przetwarzania (w sklearn metody działające w trybie online posiadają metodę partial_fit), więc w każdym kroku będziemy wykorzystywać narastające okno danych. Ponadto jako inicjalną wartość $\lambda$ i $\alpha$ będziemy podawać estymatę z poprzedniego okna.

Zdefiniujemy teraz klasę FullConjugateBayesLinReg na bazie implementacji z sklearn - sprametryzujemy ją zgodnie z naszą notacją.

Zadanie 2¶

Pewien student posiada samochód i dosyć często podróżuje z Wrocławia w Bieszczady. Jest oszczędny, więc wybiera drogę z pominięciem odcinków płatnych. Student (ani my) nie wie ile kilometrów ma ta droga, nie zna spalania swojego samochodu, a chciałby zaplanować budżet na pewną wyprawę. W ciągu ostatnich 10 podróży zanotował następujące dane:

• pieniądze wydane na paliwo [PLN] - zkładamy, że we Wrocławiu jest na stacji z pustym bakiem i dojeżdża w góry na ostatnich oparach
• średnia prędkość jazdy w danej podróży [km/h] - uwzględnia tylko czas przemieszczania się
• liczba postojów podczas danej podróży - nie zanotował jednak ile postoje trwały
• łączny czas podróży [h]

W całym rozumowaniu będziemy mieć jedno, istotne i naciągane założenie - spalanie samochodu nie zależy od prędkości jazdy i nie musi być stałe w każdej podróży.

Owy student podjął decyzję, że należy rzucić wszystko, zatankować samochód za 100 PLN i stając dwa razy na kawę jechać 90km/h aż do skończenia paliwa. O tej porze roku zachód słońca jest o godzinie 17:30. Ile maksymalnie godzin przed 17:30 powinien wyjechać żeby z 90% szansą zatrzymać się o tej porze, bądź póżniej?

Zastosuj metody poznane na zajęciach aby wyznaczyć szukaną wartość.

Zadanie 2.1 (2 pkt)¶

W komórce poniżej przygotuj dane z df do zaaplikowania następującego modelu liniowego: \begin{equation} \text{time} = \frac{\text{money}}{\text{speed}} * \text{distance_for_money_unit} + \text{num_stops} * \text{stop_time}, \end{equation}

gdzie:

• $\text{distance_for_money_unit}$ to współczynnik odpowiadający za dystans przejechany za daną kwotę
• $\text{stop_time}$ to współczynnik odpowiadający za czas spędzony na postoju

1. Utwórz zmienną 'money/speed' w zbiorze danych
2. Dokonaj standaryzacji zmiennych niezależnych, które będą wykorzystane w modelu
3. Wycentruj zmienną niezależną i zapisz jej średnią pod zmienną y_mean

Zadanie 2.2 (1pkt)¶

W komórce poniżej zainicjalizuj wybrany model Bayesowskiej Regresji Liniowej pod zmienną fuel_model i wykorzystaj go w animacji inkrementacyjnego uczenia - animacja jest w pełni przygotowana.

1.ConjugateBayesLinReg¶

2. FullConjugateBayesLinReg¶

Dane testowe na podstawie treści zadania:¶

• 100 PLN
• 90km/h
• 2 postoje

Zadanie 2.3 (1pkt)¶

Szukamy teraz wartości od, której $\text{time}$ jest większy z prawdopobieństwem co najmniej 0.9 (90% szans na dojechanie). Równoważnie wartości, od której $\text{time}$ jest mniejszy z prawdopobieństwem co najmniej 0.1. Czyli kwantylu rzędu $p=0.1$.

$$P_X((-\infty, x_p])) \geqslant p \iff P_X([x_p, \infty))) \geqslant 1-p$$

Dla rozkładów z modułu stats pakietu scipy ta funkcjonalność zaimplementowana jest pod metodą ppf.

Na podstawie rozkładu predykcyjnego dla danych testowych test_data oblicz maksymalną liczbę godzin (w sytemie dziesiętnym - 15min = 0.25h itd) zgodnie z treścią zadania i przypisz wartość do zmiennej worst_case_time.

UWAGA: pamiętaj jakie przekształcenia były stosowane na zbiorze danych!

Implementacja z dowolonym rozkładem na $\pmb{\beta}$, $\alpha$ oraz $\lambda$ z wykorzystaniem Variational Inference¶

W sytuacji gdy chcemy operować na dowolnych rozkładach prior nie jesteśmy w stanie przeprowadzić dokładnej estymacji Bayesowskiej. Ewaluacja rozkładów posterior w sposób numeryczny jest w ogólności poza zasięgiem możliwości obliczeniowych. Rozwiązaniem tego problemu jest przybliżone wnioskowanie Bayesowskie (ang. approximate Bayesian inference), gdzie miedzy innymi możemy wymienić następujące metody:

• Markov Chain Monte Carlo - poprzez błądzie losowe w oparciu o Markov Chain daje aproksymację skomplikowanych rozkładów ciągłych za pomocą techniki Monte Carlo
• Expectation Maximization - w interacyjnej procedurze odszukuje lokalne optima parametrów dla MLE i MAP gdy model zależy od pewnych nieobserwowanych zmiennych
• Variational Inference - problem maksymalizacji prawdopobieństwa brzegowego (Marginal Likelihood) jest przeniesiony na problem optymalizacji gradientowej na dolnym ograniczeniu likelihoodu - Evidence Lower BOund (ELBO)

Biblioteka pyro wspiera MCMC oraz VI, a konkretnie rozszerzoną wersję Stochastic Variational Inference, która umożliwia przetwrzania batchowe - analogicznie do Stochastic Gradient Descent.